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ABSTRACT 

As the aviation industry moves towards higher efficiency 
electrical power generation, all electric aircraft, or zero 
emissions and more quiet aircraft, fuel cells are sought as the 
technology that can deliver on these high expectations. The 
Hybrid Solid Oxide Fuel Cell system combines the fuel cell 
with a microturbine to obtain up to 70% cycle efficiency, and 
then distributes the electrical power to the loads via a power 
distribution system. The challenge is to understand the 
dynamics of this complex multi-discipline system, and design 
distributed controls that take the system through its operating 
conditions in a stable and safe manner while maintaining the 
system performance. This particular system is a power 
generation and distribution system and the fuel cell and 
microturbine model fidelity should be compatible with the 
dynamics of the power distribution system in order to allow 
proper stability and distributed controls design. A novel 
modeling approach is proposed for the fuel cell that will allow 
the fuel cell and the power system to be integrated and 
designed for stability, distributed controls, and other interface 
specifications. This investigation shows that for the fuel cell, 
the voltage characteristic should be modeled, but in addition, 
conservation equation dynamics, ion diffusion, charge transfer 
kinetics, and the electron flow inherent impedance should also 
be included. 

INTRODUCTION 

Conventional electric power generation based on coal and 
hydrocarbon fuels causes high pollution and generates power 
with relative low efficiency. In aviation and in other remote 
power applications, like in the auto industry, weight and 
volume also become important design considerations. The 
emerging fuel cell technology offers a promising solution for 
powering future aircraft with reduced emissions, higher 
efficiencies, and with lower weights and volumes. Especially, 
the emissions from the fuel cells are known to be the lowest for 
all of the fossil fuels (refs. 1 and 2). 


The hybrid Solid Oxide Fuel Cell (SOFC) power system, 
which is one of the top contenders for future aircraft power, 
incorporates a microturbine Gas Turbine Engine (GTE), a Fuel 
Cell (FC), and power processing and distribution hardware, 
shown in figure 1 . The purpose for the GTE is to improve the 
overall system efficiency by utilizing the hot air flowing out 
from the SOFC and expanding it in the turbine to generate 
additional electrical power (refs. 1 and 2). The FC and the GTE 
are tightly coupled through auxiliary components such as 
heaters, pumps, the recuperator, heat exchangers, valves, the 
reformer, humidifiers, and diffusers. In this system, the fuel 
cell, together with the heat exchanger, effectively replaces the 
combustor section in a typical gas turbine engine as shown in 
figure 2. The GTE and the FC both generate electrical energy, 
which is delivered to electrical loads by an electrical power 
processing and distribution system. The inclusion of the GTE in 
the hybrid SOFC power system improves the overall efficiency 
by approximately 20%. In a typical hybrid SOFC system, the 
FC will supply approximately 80% of the total output power 
and the GTE channel will supply 20% (refs. 2 and 3). 

The hybrid SOFC power system involves interdisciplinary 
engineering, and the approach to systems modeling for 
dynamics and controls needs to involve this type of analytical 
knowledge. Often in fluidic type systems like aircraft engines, 
the component flows can be choked. This serves to separate the 
source and load impedances of the upstream and downstream 
components, thereby, decoupling their dynamics. In electrical 
power systems, choking the current flow is not an option. 
Therefore, the source/load dynamics must be dealt with by 
designing stability margins for the interface impedances 
(ref. 4). For stability and distributed controls design, the power 
system dynamics of interest typically extend to 10’s of KHz 
(depending on the switching frequency of the power 
processors) (ref. 5). Because in this case the SOFC is the power 
source, where the power distribution system behaves as the 
electrical load (ref. 6), the source or fuel cell dynamics should 
also be modeled up to this frequency range. 
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FIGURE 1.— HYBRID SOFC POWER SYSTEM STRUCTURE 



Power turbine 

FIGURE 2.— GAS TURBINE ENGINE DIAGRAM 

This frequency range (up to 10’s of KHz) for the SOFC also 
includes conservation dynamics, diffusion, kinetics, and 
inherent impedances. Obviously, simulating the system in this 
extensive frequency range that starts from sub hertz frequencies, 
with the temperature dynamics, can not be done effectively in 
the time domain. However, these models can be converted into 
the frequency domain, which allows the applicability of 
classical control and stability design approaches. Once the 
system is designed for distributed control and stability, more 
simplified models can then be used to design and assess the 
longer time scale system operability and performance (refs. 7 
and 8). These higher fidelity models can be generated based on 
physics, which allows for more in depth understanding of how 
the process works, followed by calibration of the models based 
on test results. These models can also be generated strictly by 
testing, requiring impedance measurements in electrical systems 
or, equivalently, impedance spectroscopy measurements in 
electrochemical systems (ref. 9). 

Both the GTE compressor and turbine are directly interfacing 
with the SOFC, see figure 1. Also, the turbine through the 
electrical generator interfaces with both the SOFC and the 
electrical power management and distribution system (ref. 10). 
Therefore, the GTE model fidelity (i.e., the compressor and 
turbine) also needs to be comparable to the rest of the system by 
employing the gas dynamics through the conservation 


equations. In some cases, the flows of the GTE and the FC 
maybe choked to decouple these dynamics and in such cases the 
models can be simplified. Work to complete the GTE portion of 
this model is continuing and more detailed results may be 
presented in the future. 

This paper is organized as follows. First, details pertaining to 
the justification for high fidelity modeling of the SOFC for 
distributed controls and stability design are presented. This is 
followed by presenting a physics based model of the SOFC, 
suitable for power systems controls and stability design. Finally, 
concluding remarks are presented. 

FUEL CELL MODELING IN REMOTE ELECTRICAL 
POWER DISTRIBUTION SYSTEMS 

In general, (without distinction of electrical, mechanical, or 
electrochemical systems) when an upstream component 
provides energy processed or utilized by a downstream 
component, the two components connected together form a 
“source-load” configuration as shown in figure 3 for a power 
system (ref. 4). In this figure, G s (s ) and G L (s ) stand for the 
corresponding input to output transfer functions, Z 0 (s) and Z L (s) 
stand for the source output impedance and the load input 
impedance respectively, R L represents the load resistance, and 
Zeq, Zp are the source series impedance, and the load power 
stage impedance respectively, which can also include filters. If 
these impedances are sufficiently separated as a function of 
frequency (>6dB or twice the magnitudes separation), then the 
overall transfer function formed by connecting these 
components can be approximated by the product of the 
individual component transfer functions. However, if this 
condition is not met, then the overall transfer function is 
described by the following equation: 


G t ( s ) = 


G S {s)G L {s) 


1 + 


Zp(s) 

Z L (s) 


( 1 ) 
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FIGURE 3.— FUEL CELL AND POWER DISTRIBUTION 
EQUIVALENT INTERFACE NETWORK MODEL 


As seen, the denominator of equation (1) forms a characteristic 
equation, whose Nyquist stability criteria need to be defined and 
designed for. Inspection of equation (1) shows that if the 
magnitude of Z 0 (s)/Z L (s) becomes 1 at some point in frequency, 
while its phase becomes ±180°, the overall transfer function will 
be infinite, and thus completely unstable. There is a range of 
magnitudes and phases around this point where the overall 
transfer function of equation (1) will exhibit oscillatory 
behavior, which is also undesirable. The output impedance Z 0 (s) 
depends on the hardware designs and physics, of the fuel cell 
components in this case, as well as any feedback control 
designs. 

Typical electrical power systems do not meet this condition of 
sufficient separation in impedances, which makes it necessary to 
design for adequate stability margins. Adequate stability 
margins are typically 60° phase margin and 6 to 10 dB gain 
margin. Designing adequate stability margins is vitally 
important in remote power system applications, which are not 
afforded the stability of large interconnected networks like the 
terrestrial electrical power network. It is also important to 
design for adequate stability margins in systems with increased 
constant power or negative impedance loading. This is also true 
for remote electrical power systems where the load closely 
matches that of the source generating capacity, like future 
aircraft electrical power systems. 

Figure 4 shows a generic Nyquist stability diagram for this 
characteristic equation, and a designated exclusion zone 
designed such that the interface is specified with absolute and 
conditional stability margins of >60° phase and 6 dB gain. In 
other words, the magnitude of the Z 0 (s)/Z L (s) ratio is not 
permitted to become greater than 0.5 while the phase of this 
ratio resides between 120 and 240°. Given the output impedance 
of the FC, the input impedance of the power distribution system 
(or the inverter in this case) interfacing the FC, will need to be 
designed so that the Nyquist stability criteria is satisfied. 
Figure 5 shows the power system stability interfaces for this 
simplified hybrid SOFC system diagram. Instabilities like 
electrical instabilities in this case (depending on their severity); 
can propagate through the system to also involve mechanical 
components, which can cause catastrophic system failures. 

The following is a description of general requirements for 
distributed controls and stability design that pertain to interface 
impedances and closed loop controls design. For interface 
impedances, in addition to the Nyquist criterion, all frequencies 
should be damped. For closed loop distributed controls, the 



Real axis 


FIGURE 4.— NYQUIST STABILITY OF Z 0 (s)Z L (s) FOR 
ABSOLUTE AND CONDITIONAL STABILITY OF > 60° 
PHASE MARGIN AND > 6 DB GAIN MARGIN 



FIGURE 5.— TYPICAL HYBRID SOFC POWER SYSTEM 
STABILITY INTERFACES (S') 

frequencies of the various control loops should be separated 
when feasible, and adequate damping, phase and gain margins, 
and disturbance attenuation should be designed. For power 
stages, in addition to distributed controls, design to satisfy the 
Middlebrook (ref. 11) stability criterion and audio susceptibility. 
The design for all these specifications will necessitate 
developing high fidelity models of the type discussed here. For 
completeness, there is also stability implications associated with 
any device that produces thermal energy; for instance, the 
energy produced by the GTE and the FC could be susceptible to 
thermo-acoustic and the Helmholtz instabilities, but these can be 
addressed separately. 

SOFC MODELING FOR SYSTEM DYNAMIC 
INTERACTIONS, CONTROLS AND STABILITY 

In the previous section, the justification for the need of high 
fidelity dynamic modeling for the FC subsystem was covered. It 
was established that the FC dynamics need to extend into 10’s 
of KHz to make it comparable to the dynamics of interest in a 
typical electrical power distribution systems for control and 
stability analysis. This frequency bandwidth for the FC includes 
the following areas for modeling: The FC voltage relation, 
conservation equations modeling for species, momentum, and 
energy that covers sub hertz up to 100’s of Hz, ion diffusion that 
cover 100’s of Hz, charge transfer kinetics that can extend to 
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10’s of KHz, and finally the inherent impedance of the fuel cell 
relative to the flow of electrons. 

The conservation equations govern the mass transport 
phenomena into the FC and the temperature and flow rates of 
the species. The gas diffusion governs the transport of ions 
through the porous electrode and electrolyte media to the 
reaction site. This reaction site is called “triple point boundary” 
( tpb ), which consists of the interface between the electrode, the 
electrolyte and the gasses that react. The charge transfer kinetics 
governs the electrochemistry reactions at the tpb associated with 
the charge transfer and charge balance. 


Fuel Cell Voltage 

In the ideal gas law, the electromotive force (EMF) or 
reversible open circuit voltage is related to the Gibbs free 
energy (refs. 1 and 2) release as 


E = - 


A G 
nF 


(3) 


where AG is Gibbs free energy, F is Faraday’s constant [96487 
C mol -1 ], and n is the number of electrons participating in the 
reaction, which is 2 for the hydrogen/oxygen reaction. This open 
circuit voltage is related to the gas partial pressures and 
temperature through the following equation: 


E = E° 


RT t 

+ In 

2 F 


Pk 


Ph 2 o 


(4) 


where E° is the voltage associated with the reaction free energy, 
R is the universal gas constant [82.05E-5 J mol -1 K _1 for 
pressure in atm], T is the temperature, pj are the partial 
pressures of the associated gases at the tpb where the reactions 
take place, as shown in figure 6 (ref. 18 revised). In 
equation (4), it is assumed that the vapor pressure of the steam 
at the temperature concerned ( p°H 20 ) is equal to the standard 
pressure (p°). If this is not the case, then the partial pressures of 
hydrogen and oxygen are divided by p° , and the partial 
pressure of steam is divided by p^o • The actual cell voltage is 
reduced by irreversible losses associated with activation loss, 
ohmic loss, and concentration loss respectively (refs. 1 and 2) as 



f • \ 

z 


f • A 

V = E-A In 

- iR in -B\n 

i-i- 


\ l o ) 

l n) 


where i is the operating current of the cell, i 0 is the exchange 
current, // is the limiting current at which the fuel is used up at a 
rate equal to its maximum supply rate, R in is the inherent 
resistance of the fuel cell, and A and B are constant coefficients 
determined experimentally. There is a certain I-V characteristic 
associated with every fuel cell, which can be constructed from 
experimental data, as shown in figure 7. 

Equation (5) is a steady state equation that contains no source 
impedance dynamics. Ultimately, the objective will be to 
incorporate the appropriate dynamics so that this power source 
can be expressed as an equivalent circuit of voltage and 
impedance as shown in figure 3 and equation (1). Modeling the 
appropriate dynamics of this voltage source will allow the 
application of typical power system distributed controls and 
stability design (ref. 4). 
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FIGURE 6.— PRINCIPLE OF THE SOLID OXIDE FUEL CELL 



FIGURE 7.— TYPICAL FUEL CELL l-V CHARACTERISTIC 


Conservation Equations Modeling 

Conservation equations modeling involves mass conservation, 
energy conservation, and momentum conservation. These 
equations are expressed as one dimensional, lumped volume 
models. The conservation equations deal with the pressures, 
temperature and flows in and out of the fuel cell, shown in 
figure 6. The frequencies involved in these equations start from 
sub hertz up to approximately 100 Hz. 


dt 


RT 


f \ 

wj n - wf - wf ut 
V / 


Mass (6) 


dT_ 

dt 


1 

mC 


m_n_ 

Q — W + Wri h ri ~ W pi h pi 
V i = 1 i = 1 J 


Energy 


(7) 


d w 
dt 



Momentum (8) 


For the mass equation (13), p^ is the bulk pressure due to that 

particular specie i in the flow channel of the fuel cell. V c is the 
channel volume, and it could be different for the anode and the 
cathode. w t is the flow rate of the specie i in [mol s” 1 ] of the 

incoming flow ( wf 1 ) , the flow that diffuses and reacts (w- ) and 
the outgoing flow {w ° ut ) , respectively. The reacting species 
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considered with the SOFC are //?, H 2 O f and 0 2 . For the 
development in this subsection, ion diffusion is considered as a 
steady state process. Then the reactive flow in equation (6) can 
be represented as 


Ngl 

nF 


(9) 


where N 0 is the number of cells in the stack in series, / is the 
stack current, with n = 2 for hydrogen and water and n = 4 for 
oxygen (in consideration that for the same number of electrons 
involved 1/2 mole of oxygen reacts for every mole of 
hydrogen), figure 6. Later, ion diffusion will be treated for non- 
steady state conditions. It could be considered that the molar 
flow of any gas through the valve is proportional to the partial 
pressure inside the channel (ref. 12) as 


Wi=K iP b (10) 

where K t is the valve molar constant for the particular specie. 
Equation (10) can then replace w° ut in equation (6). 

For the energy equation (ref. 13), additional species could be 
considered of atmospheric air C0 2 for the water gas shift 
reaction. This energy equation uses steady state cell 
electrochemistry calculations based on a Microsoft Excel model 
developed by Pacific North National Laboratory (ref. 14) 
(PNNL) that was converted here into MATLAB®. In this PNNL 
model, temperature of the cell is assumed to be constant, 
representing the cell in a controlled temperature environment. 
This has been changed here, as it is not the case with a fuel cell 
stack, which adds heat loss to the model and a non fixed exit 
temperature. This change was incorporated based on 
equation (7), which adds temperature dynamics to the model. 
For the energy equation, T signifies the cell or the stack 
temperature (i.e., the stack by dividing the equation by the 
number of cells in the stack). The mass is m [kg], C is the 
specific heat capacity [J kg -1 K 1 ], which accounts for the heat 
capacities and the mass of the individual fuel cell materials, see 
table 1 for typical fuel cell materials (ref. 15). Subscript r 
signifies the reactants while subscript p signifies the products. Q 
is the heat transfer rate [J s _1 or watts], which is the difference 
between the fuel cell open circuit voltage equation (3) and the 
operating voltage times the current in figure 7, W is the 

electrical work transfer rate [watts], and h is the total enthalpy 
per molar basis [J mol -1 ], which is the sum of the enthalpy of 
formations and the change in enthalpy from atmospheric 
conditions. It is assumed that the only heat loss is through the 
exit stream, because the fuel cell stack is well insulated so that 
the radiating heat loss is considered to be negligible. 


TABLE 1.— SOFC SINGLE CELL DATA 



Spec, heat 
cap. (J/kg*K) 

Mass Density 
(kg/m 3 ) 

Thickness 

(m) 

Anode 

466 

4200 

6x10“ 

Cathode 

520 

6350 

5x10 s 

Electrolyte 

460 

6010 

1x10 s 

Interconnect 

800 

7700 

5x1 0" 4 
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FIGURE 8.— FUEL CELL STACK SIMULATION FOR A STEP IN 
LOAD CURRENT (TEMP, VOLTAGE and PARTIAL PRESS) 


For the momentum equation, A is the flow channel cross section 
area [m\ l is the length of the channel \m\, g is the gravitational 
acceleration [9.8m s“ 2 ], and p T is the total pressure [atm]. 

The conservation equations (6 to 8) tie to the voltage 
equation (5) by the temperature and partial pressures in 
equation (4). If ion diffusion is assumed to be a steady-state 

process, then the partial pressures at the bulk flow p 1 - can be 

substituted for pj in equation (4). If ion diffusion dynamics are 

taken into account, as will be covered in the next subsection, 
then the partial pressures at the tpb are computed from the bulk 
partial pressures in equation (6). Figure 8 shows a simulation of 
a fuel cell stack involving the conservation equation dynamics 
and the fuel cell voltage relation for a step in load current. 
Momentum is not simulated presently due to the specifics of 
stack design associated with geometric properties. 

Ion Diffusion Modeling 

Ion diffusion deals with the diffusion dynamics of the gas 
species starting from the channel flow stream entering the 
electrode and electrolyte porous material to the tpb where the 


electrochemical reaction occurs, shown in figure 6. In 
equation (5), there is a voltage loss attributed to the 
concentration of species. In addition, there are diffusion 
dynamics associated with species concentration and partial 
pressures. These dynamics do not produce a loss per se, but 
rather they contribute to the voltage dynamics as a reactive 
impedance of the fuel cell materials to the flow of ions. The 
diffusion dynamics are expected to cover a frequency range 
extending from 100’s of Hz up to a few kHz. 

The diffusion equation, which in this case describes the 
species or mass transport from the flow stream to the tpb 
reaction site, can be written as 



d 2 C 
dx 2 


( 11 ) 


where C is the species concentration, D is the effective diffusion 
coefficient, and x is the diffusion depth as shown in figure 6. For 
more detailed calculations of the diffusion coefficient based on 
the porosity and tortuocity of the fuel cell materials, see (refs. 16 
and 17). One way to solve equation (1 1) is to perform a Laplace 
transformation to convert the partial differential equation into an 
ordinary differential equation (18). 

-iow-0 (.2) 



with the boundary conditions 

^^ = -Tj o (V)| x=0 , CHs) = C(s) \ X=L (13) 
dx D 

where J Q is the diffusion flux [mol s _1 cm -2 ]. Solving 
equation (12) in terms of C(s,x ) results in an expression of 
exponential functions, which can be represented in transfer 
function form by applying Taylor series expansion (ref. 18). 
Also applying the ideal gas law, this results in the following 
equations: 

Pi 00 = GjpJj (s) + Gpppf (. s ) (|4 . 

Jl’(s) = Gjjj;(s) + G PJ p'’(s) ( 

where 


jrx G N 

(JJP = K JP — 3 

Gd 


Gpj = K P j 


Gd 


°pp = -77- > 

Gd 


Kjp = - 


RTL 

~AD 


r Ll 4.1 
Cjm = S + 1 , 

6D 


Gjj = Gpp , 
AL 


K PJ =- 


RT 


L 2 


Gd - 1 + + 


Z 4 


2D 24 D 2 


with the partial pressures p at the tpb and the bulk flow 
represented by the superscripts t and b respectively, f and J 
signify the diffusing flow rate at the bulk flow boundary and the 
reactive flow rate at the tpb respectively in [mol s -1 ]. The 
diffusion cross section area is A [m 2 ], and L is the thickness of 
the diffusion layer [m], which can be different for the anode and 
the cathode. In figure 6, it can be seen that oxygen in addition 


diffuses through the electrolyte layer. However, because the 
electrolyte layer is ionically highly conductive, its associated 
dynamics could be much faster, and it is assumed that these 
dynamics can be neglected. Inspection of the transfer functions 
in equation (14) shows that that their frequencies are as follows 


0) p =2yf6-^, 


©Z 1 = 


6D_ 

I? 


where ($ p , co z i are the frequency of the double pole location and 
the zero, respectively [rad s -1 ]. The zero of G PJ is located at zero 
radians. Based on this, the diffusion response time will be 
proportional to the D/L 2 ratio, which means that the response is 
inversely proportional to the square of the diffusion layer 
thickness and proportional to the effective diffusion constant. 

Similarly (as in the previous subsection), the assumption here 
is that the charge transfer kinetics is a steady state process in 
order to allow coupling of the diffusion process to the 
conservation equations. This assumption implies that as long as 
the reaction at the tpb is not concentration limited, the reaction 
rate [moles s” 1 ] will adjust instantaneously to meet a change in 
electrical current demand. Then the reactive flow at the tpb can 
be expressed as 

N I 

J\=^f (15) 

nF 

This relaxes the previous assumption of a steady state diffusion 
process, reflected in equation (9), and equation (9) now becomes 


W[=p (16) 

which effectively ties ion diffusion to the dynamics of the 
conservation equations. Figure 9 shows a cell simulation of the 
mass flow rates of hydrogen and oxygen at the bulk flow 
channels and their respective partial pressures at the tpb, due to 
a step in load current at 0.05 s. Also, figure 10 shows the same 
variables plotted for a step change in the bulk flow partial 
pressures. In these simulations, diffusion has not been tied to the 
conservation equations dynamics. In the more detail plots of 
figures 9 and 10, it can be seen that the steady-state molar flow 
rate of H 2 is twice that of 0 2 , which satisfies the chemical 
reactions of these species as shown in figure 6. 

Charge Transfer Kinetics Modeling 

Charge transfer kinetics, or what is also commonly referred to 
as charge balance, deals with the actual reactions that take place 
at the FC tpb' s, see figure 6. This area of research is still in its 
relative early stages in terms of understanding the reaction 
mechanisms that take place at these tpb' s and how to model 
these reactions (refs. 19 and 20). Models can also be generated 
or optimized by performing impedance spectroscopy testing, but 
testing alone would make it difficult to separate the cell 
dynamics due to diffusion from those due to kinetics. Also, 
models generated by impedance spectroscopy alone will need to 
be revaluated at various operating conditions, which can 
complicate the controls development. 
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FIGURE 9.— CELL SIMULATION OF MASS FLOW RATE AT 
THE BULK CHANNEL AND PARTIAL PRESSURES AT tpb 
DUE TO A STEP IN LOAD CURRENT 
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FIGURE 10.— CELL SIMULATION OF MASS FLOW RATE 
AT THE BULK CHANNEL AND PARTIAL PRESSURES 
AT tpb DUE TO A STEP IN PARTIAL PRESSURE AT THE 
BULK FLOW 


First some chemical background; the cell electrochemical 
model describes a set of possible chemical and electrochemical 
reactions that are derived from basic surface electrochemistry. 
These reactions are of the following form 

k] k_\ 

aA + bB < > cC 

r , (17) 

^A = k_ x a[c] C -ha[A] a [B] b 
dt 

where A, B , and C are gas phase species, surface species 
absorption sites, or electrons; & 7 and &_ 7 are reaction rate 
constants for the forward and the backward reactions, 
respectively. The chemical and electrochemical equations canbe 
formulated as mass and charge balance. The concentration or 
mass balance concerning the time dependant surface species A is 
given by the corresponding differential equation in 
equation (17). Both the anodic and the cathode reactions need to 
be considered separately. 


Modeling of the Anodic Reactions. — The predominant 
electrochemical reactions (ref. 19) concerning the SOFC at the 
anodic tpb are: 

k\ k_\ 

H2 (g) + 2 ad <■ 2 > 2 H a d 

H 2 0 (g) + ad <- h 2 0 ad 

Oad + Had* — — HO ad + ad (18) 


k/\ k_A 

H 2 O ad +0 ad * ^ ^ OH a d 
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HlO ad + ad 


* OH ad +H ad 


Tt _ N o i K 

J i „Z7 


(26) 


k^k-s 


0$ + ad <r- 


hJt-6 


> O ad + To + 2 <?M 


where g denotes the gas phase and ad denotes the absorbed 
phase. k t are the reaction rate constants for the different reaction 
steps, and O x denotes the oxygen interstitial and F an oxygen 
vacancy. Based on that, the time dependant concentration or 
fraction of surface coverage 0 f is (ref. 19) 

dO 

= 2*1 p ‘ H2 (1 - 0) 2 - 2 *_ 10 ^ - * 307/00 (19) 

+ k-?,QoH (1 - 0) + *5 9 / 7 2 0 (1 — 0) — *-5 00 / 70/7 


— = *3 0(90// ~ *-300// (1 “ 0) + 2*40// 2 O0O 
- 2k_ 4 Qj )H +*50// 2 o(1 — 0) — *-50O//0// 

= k 2P‘ H2 o 0 ~ 0) ~ *- 2 6 // 2 o ~ * 4 9// 2 o9o 

+ *-40o// -*50// 2 O 0-0) + *-500/7 0/7 


d^o 

dt 


+ £-40^ + /£- 


( 22 ) 


These changes essentially tie the conservation equations and ion 
diffusion to the charge transfer kinetics. The overpotential in 
equation (25) can be calculated as 

(27) 

where E 0 is the fuel cell open circuit voltage and F is the actual 
operating voltage of the fuel cell, equation (5). The multitude of 
reaction rate constants makes it difficult to produce meaningful 
simulation results, without the benefit of optimizing the model 
based on test data. 

For the moment, this model development assumes that the ion 
concentration of 0$ in the YSZ electrolyte is constant, thereby, 

ignoring the kinetics at the cathode tpb. However, similar to the 
development of the kinetics at the anode tpb, the kinetics at the 
cathode tpb also need to be considered. 

Modeling the Cathode Reactions. — The electrochemical 
reaction of oxygen reduction at the cathode tpb is described in 
(ref. 20). The overall reaction is considered to consist of several 
elementary reactions, mainly those associated with oxygen 
absorption, desorption, diffusion, and electronation. It is 
assumed that this process is dominated by the following 
reactions; (ref. 20) 


kad -fades 

02(g) +2 s <— — — > 2 O a d (28) 


Equations (19 to 22) are nonlinear and would need to be 
linearized by an appropriate method for frequency domain 
analysis. The total fraction of surface coverage, 0 is 
0 = 0#+ 0 oh + ®h 2 o + ®o- The reaction rate constants of and k 2 
in k\p l H and k 2 P t H ^ Q are not constant as the partial pressures 

are modified. The interface reaction is 

OS + (1 - 0) < k,Jc -' > 0 O + K + 2 ejfi (23) 


O ad + V” + 2e- 4- - ,C —> 0$ +s (29) 

where s is the concentration of vacant surface sites, similar to 
(1-0) from before and k ad , k des signify the oxygen absorption and 
desorption reaction rates. If the model is normalized by 
introducing N 0 , the maximum number of surface sites per unit 
area which can be occupied by O a d , then the time dependent 
concentration or fraction of surface coverage 0 a d •> ( re T 20), is 


Electrical current is related to the rate of electron production 
crossing a given cross section area per second, which is related 
to the molar flow rate of the reaction as 

._dQe- _ 1 dVS _ I dQ 0 _ 1 _ dOS _ 1 +1-9) (24) 

dt 2 dt 2 dt 2 dt 2 dt 

where q is the electrical charge [Coulombs]. Current flows in 
the direction of e production as i-if- i r • Applying the Buttler- 
Volmer equation and interface reaction gives 

— 2(1 — P)(— — — r|) -20 (— ti) 

i K =k 0 {e PA RT v [ 0 g]( l-Q)-e RT 9 Q 0 [V 0 ']} (25) 

where k 0 is a reaction constant, p is a charge transfer coefficient, 
and r| is the overpotential, see figure 7. In equations (19) and 
(21) the partial pressures at the tpb are substituted in the place of 
the partial pressures at the bulk flow, used in (ref. 18). With the 
inclusion of the kinetics modeling of equations (19 to 22), and 
the calculation of the charge balance in equation (25), 
equation (15) in the previous section becomes 


= 2k adPo No (\ - Qad) 2 - 2 k des No& ad 
-i c l(2N 0 FA c ) 


i c = 2 N 0 FA C [*_ lc [OS ](l -9 lu l) ~ he [vj]Q ad 


(31) 


where the reaction rate constants k lc and k. lc are 

k-\ c = k-\cOe 


F 

-2p T| 

k\ c = hcoe RT 


-2(l-p)+ti 

RT 


and k lc0 and k. lc0 are potential independent reaction rates that 
depend on temperature. Alternatively, equation (25) can also be 
expressed exactly as equation (31) where it is a function of its 
triple point boundary area with the reaction rates of k 6 

substituted for k lc . The time dependent concentration of Oq , 

could be expressed as 


NAS A/TM— 2006-2 14104 



(32) 


10 


Bode diagram 


— k\c®ad k-lc^ox 
Oq ~ Pox^ox 


(33) 


where p ox is the maximum density of reduced oxygen per unit 
area [moles-m“ 3 ]. Equations (32) and (33) are an extension of 
that presented in (ref. 20) and may need to be further 
investigated. Equation (33) can be substituted into (25) to 
calculate the anodic current. The anodic and the cathode 
currents in equations (25) and (31) will need to be balanced, 
which will call for an iterative procedure. The overpotential rj 
can be different for the anode and the cathode based on the 
voltage equations for their respective tpF s, and the overall 
voltage would be the difference of the anodic and the cathode 
voltages as 

V = Va-Vc (34) 

The anodic and the cathodic impedances can be computed 
separately based on the models presented above and the reaction 
rate constants (refs. 19 and 20) at some operating point by 
applying a sinusoidal perturbation as 


Z 


di 


dp 


- oper 


(35) 


Figure 1 1 shows the cathodic impedance plot of equations (30), 
(31), and (35) for an operating overpotential of 0.1 V, for certain 
temperature and pressure, and without combining the anodic and 
cathodic systems. The anodic impedance has about the same 
shape as the cathodic impedance. But the anodic impedance is 
more difficult to determine without the aid of test data due to the 
multitude of reaction rate constants. 

The FC reactions are still an area of intensive research. The 
reactions presented here may not describe too accurately what 
happens at these tpb' s. However, in the case of the anodic 
kinetics, which in this case is described by a fourth order model, 
it is expected that the order of the model is high enough to 
accurately model the impedance data obtained by testing. 


Fuel Cell Impedance 

Inherent impedance of the FC is complex, and there are many 
references that devote discussions in this area 
(refs. 9, 18, 19, and 20). Figure 12 shows an equivalent circuit 
of the FC impedance. In this circuit, Cdi a and Q/ c represent the 
charge double layer capacitance between the anode, cathode and 
the electrolyte. The anode and cathode diffusion impedances, 
which are not explicitly shown here, are represented by Z da and 
Z dc • The combination of R t , R p , and Cj form the kinetic 
impedances for the anode and cathode. The charge transfer 
resistance is represented by R t and R p stands for the polarization 
resistance, and C; is an equivalent capacitance of the kinetic 
impedance. The electrolyte resistance is represented by R e . In 
this circuit the resistances of the anode and the cathode 
electrodes together with the interconnect resistances are not 
shown. But if these resistances are significant compared to the 
rest of the impedance circuit, they can be added to the respective 
ends of the circuit in figure 12. For complete FC impedance for 
the purposes of control and stability analysis, the equivalent 



Frequency, rad/sec 

FIGURE 11.— BODE PLOT OF THE CATHODIC IMPEDANCE 


Cdlc Cdla 



Clc Cla 


FIGURE 12.— EQUIVALENT FUEL CELL CIRCUIT 


impedances due to the conservation equations (6 to 8), plus 
those of any feedback controls would need to be included. This 
work is still in progress. 

These impedances can be computed by linearizing the models 
in MATLAB®, computing the state space matrices, and then 
plotting the response in a Bode plot. The equivalent circuit can 
then be constructed based on the Bode plot. For the Bode plot of 
figure 11, the low frequency intercept is the resistance R pc +R tc 
in the circuit of figure 12. The low frequency comer frequency 
is cd zc = l/[(R p +R t )C!] in radians/sec. The higher frequency 
intercept is equal to R tc , and the upper comer frequency is 
cDp C =l /(RfCi). Based on this R pc ~ 3 Cl, R tc ~ 0.2 Cl and C lc ~ 
5.6x1 0~ 4 F. As a result, the cathodic impedance would be 

Zkc(s) = [Rpc + + (36) 

\S / (Op C + 1 ) 

The output impedance of the circuit, which is important for 
stability design of the fuel cell power system as was discussed 
earlier, can be computed by applying a sinusoidal small 
magnitude current sweep at the right hand side of the circuit of 
figure 12 and computing the voltage drop due to this excitation 
as shown in figure 13 as 




dV 

~di 


-ioper 


(37) 
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-Zdc 


Cdlc 

-Hh 


RtC ^aaA-, ~^ Sr ~ HVVVn Rta 

Uvvcpj Lq^v-J 


Cdla 

He- 


Zda 


Clc 


Cla 


J 


t ^ 

Y 0/ 


FIGURE 13.— CIRCUIT FOR COMPUTING OUTPUT 
IMPEDANCE 


To actually compute the FC output impedance the conservation 
equation equivalent impedance will also need to be included, 
which would also depend on any regulating control loops for the 
fuel cell. 

CONCLUSION 

Modeling the Solid Oxide Fuel Cell to compute its impedance 
is important in the area of electrochemistry in order to 
understand and improve the fuel cell performance. As covered 
in this paper, this area of modeling is also important in order to 
understand how to design to interface the fuel cell with remote 
power systems for applications such as those found in aviation, 
terrestrial vehicles, and the ship power industry. For this 
purpose, generic models have been developed to cover the fuel 
cell conservation equations, ion diffusion, charge transfer 
kinetics, and inherent impedance modeling. So far, most of the 
models covered in this paper have been presented in the time 
domain. 

To develop the fuel cell impedance for control and stability, 
these models will need to be incorporated in the frequency 
domain with equivalent circuit representations. As such, the 
work to develop the overall frequency domain model of the fuel 
cell impedance is continuing, as well as calibrating the models 
with actual fuel cell test data. It may turn out through testing 
that more than one frequency domain model is needed to cover 
the fuel cell operating range. Also, for a complete model of the 
hybrid FC power system the gas turbine model, the power 
distribution, and the various auxiliary component models will 
need to be developed. This work is still in progress in various 
degrees of completion. 

REFERENCES 

1. National Energy Technology Laboratory, “Fuel Cell 
Handbook 6th Edition, “ Morgantown, WV, 2002. 

2. Larminie J.; Dicks A.; “Fuel Cell Systems Explained,” John 
Wiley and Sons, Inc., New York, NY, 2003. 

3. Schonewald R.; “Gas Turbine Technology for SOFC-GT 
Systems,” International Colloquium in Environmentally 
Preferred Advance Power Generation, ICEPAG, 2004. 

4. Lee F.C.; Boroyevich D.; “Modeling and Control Design of 
DC/DC Converters,” Short Course in Power Electronics, 
Virginia Polytechnic and State University, Blacksburg, VA. 


5. Ridley R.B.; “A new Small-Signal Model for Current-Mode 
Control,” Proceedings of Power Conversion and Intelligent 
Motion Control Conference, 1989. 

6. Wang C.; Lu B.; Qiu Y.; Huang B.; Xu M.; Lee F.C.; 
Kopasakis G.; “Control of the Hybrid Solid Oxide Fuel Cell 
Power System,” CEPS Conference, Blacksburg, VA, 2005. 

7. Huang B.; Qiu Y.; Lu B.; Wang C.; Xu M.; Lee F.C.; 
Kopasakis G.; “Modeling and Simulation of the Start-up 
Process of the Turbo Subsystem in Hybrid SOFC Power 
Systems,” CEPS Conference, Blacksburg, VA, 2005. 

8. Qiu Y.; Lu B.; Wang C.; Huang B.; Xu M.; Lee F.C.; 
Kopasakis G.; “Modeling and Control of Turbine/Generator 
in the Solid Oxide Fuel Cell Hybrid Power Systems,” CEPS 
Conference, Blacksburg, VA, 2005. 

9. Macdonald, J.R.; “Impedance Spectroscopy,” John Wiley & 
Sons, Inc., New York, NY, 1987. 

10. Xu M.; Wang C.; Qiu Y.; Lu B.; Lee F.C.; Kopasakis G.; 
“Controls and Simulation for Hybrid Solid Oxide Fuel Cell 
Power Systems,” Applied Power Electronics Conference, 
Dallas, Texas, 2006 (pending). 

11. Middlebrook, R.D.; “Topics in Multi-Loop Regulators and 
Current-Mode Programming,” Proceedings of IEEE Power 
Electronics Specialist Conference, June, 1985. 

12. Padulles J.; Ault G.W.; “An Integrated SOFC plant 
Dynamic Model for Power Systems Simulation,” Journal of 
Power Sources, 86, 495-500, 2000. 

13. Bejan A.; “Advanced Engineering Thermodynamics,” John 
Wiley & Sons, Canada, 1997. 

14. Chick L.A.; Williford R.E.; Stevenson J.W.; Windisch Jr. 
C.F.; Simner S.P.; “Experimentally-Calibrated Spreadsheet- 
Based SOFC Unit-Cell Performance Model,” Proceedings 
of Fuel Cell Seminar, Palm Springs, CA, 2002. 

15. Khaleel M.A.; Lin Z.; Singh P.; Surdoval W.; Collin D.; A 
Finite Element Analysis Modeling Tool for Solid Oxide 
Fuel Cel Development,” Journal of Power Sources, 130, 
136-148, 2004. 

16. Reid R.C.; Prausnitz J.M.; Poling B.E.; “The Properties of 
Gases and Liquids,” 4th Ed., McGraw-Hill, New York, 
1987. 

17. Suwanwarangkul R; et al.; “Performance comparison of 
Fick’s, Dusty-Gas and Stefan-Maxwell Models to Predict 
the Concentration Overpotential of a SOFC Anode,” 
Journal of Power Sources, 122, 9-18, 2003. 

18. Qi Y.; Huang B.; Chuang K. T.; “Dynamic Modeling of 
Solid Oxide Fuel Cell: The Effect of Diffusion and Inherent 
Impedance,” Journal of Power Sources, 150, 32-47, 2005. 

19. Bieberle A.; “The Electrochemistry of Solid Oxyde Fuel 
Cell Anodes: Experiments, Modeling, and Simulations,” 
Ph.D. Dissertation, Swiss Federal Institute of Technology, 
2000. 

20. Mitterdorfer A.; “Identification of the Oxygen Reduction at 
Cathodes of Solid Oxide Fuel Cells,” Ph.D. Dissertation, 
Swiss Federal Institute of Technology, 1997. 


NAS A/TM— 2006-2 14104 


10 


REPORT DOCUMENTATION PAGE 

Form Approved 
OMB No. 0704-0188 

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson 
Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503. 

1. AGENCY USE ONLY ( Leave blank) 

2. REPORT DATE 

December 2006 

3. REPORT TYPE AND DATES COVERED 

Technical Memorandum 

4. TITLE AND SUBTITLE 

5. FUNDING NUMBERS 


A Theoretical Solid Oxide Fuel Cell Model for System Controls 
and Stability Design 


6. AUTHOR(S) 

George Kopasakis, Thomas Brinson, Sydni Credle, and Ming Xu 


WBS 599489.02.07.03 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
John H. Glenn Research Center at Lewis Field 
Cleveland, Ohio 44135-3191 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 


E- 15439 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 

NASA TM - 2006-2 14104 
GT2006-91247 


11. SUPPLEMENTARY NOTES 

Prepared for the Turbo Expo 2006 sponsored by the American Society of Mechanical Engineers, Barcelona, Spain, 
May 8-11, 2006. George Kopasakis, NASA Glenn Research Center; Thomas Brinson and Sydnia Credle, Florida 
A & M University, 1500 Wahnish Way, Tallahassee, Florida, 32307; Ming Xu, Virginia Polytechnic and State 
University, 118 North Main Street, Blacksburg, Virginia 24061. Responsible person, George Kopasakis, organization 
code RIC, 216-433-5327. 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified - Unlimited 

Subject Categories: 07, 20, 25, and 34 

Available electronically at http://gltrs.grc.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 301-621-0390. 


12b. DISTRIBUTION CODE 


1 3. ABSTRACT (Maximum 200 words) 

As the aviation industry moves towards higher efficiency electrical power generation, all electric aircraft, or zero 
emissions and more quiet aircraft, fuel cells are sought as the technology that can deliver on these high expectations. 

The Hybrid Solid Oxide Fuel Cell system combines the fuel cell with a microturbine to obtain up to 70 percent cycle 
efficiency, and then distributes the electrical power to the loads via a power distribution system. The challenge is to 
understand the dynamics of this complex multi-discipline system, and design distributed controls that take the system 
through its operating conditions in a stable and safe manner while maintaining the system performance. This particular 
system is a power generation and distribution system and the fuel cell and microturbine model fidelity should be 
compatible with the dynamics of the power distribution system in order to allow proper stability and distributed controls 
design. A novel modeling approach is proposed for the fuel cell that will allow the fuel cell and the power system to be 
integrated and designed for stability, distributed controls, and other interface specifications. This investigation shows 
that for the fuel cell, the voltage characteristic should be modeled, but in addition, conservation equation dynamics, ion 
diffusion, charge transfer kinetics, and the electron flow inherent impedance should also be included. 


14. SUBJECT TERMS 

Fuel cells; Modeling; Dynamics; Systems stability and control 


15. NUMBER OF PAGES 

16 


16. PRICE CODE 


17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 


18. SECURITY CLASSIFICATION 
OF THIS PAGE 

Unclassified 


19. SECURITY CLASSIFICATION 
OF ABSTRACT 

Unclassified 


20. LIMITATION OF ABSTRACT 


NSN 7540-01-280-5500 


Standard Form 298 (Rev. 2-89) 
Prescribed by ANSI Std. Z39-18 
298-102 



